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Abstract 

In the absence of losses the phase of a Bose-Einstein condensate undergoes collapses 
and revivals in time due to elastic atomic interactions. As experiments necessarily involve 
inelastic collisions, we develop a model to describe the phase dynamics of the condensates 
in presence of collisional losses. We find that a few inelastic processes are sufficient to 
damp the revivals of the phase. For this reason the observability of phase revivals for 
present experimental conditions is limited to condensates with a few hundreds of atoms. 
PACS numbers: 03.75.Fi, 42.50.Gy, 05.30.Jp 

1 Introduction 

Since the recent experimental observations of Bose-Einstein condensation in dilute atomic gases 
[|TJ, ||, ||, much interest has been raised about the characteristic features of the condensate 
[|J, and about its coherence properties in particular. Considerable attention has been devoted 
to the matter of the relative phase between two Bose-Einstein Condensates (BECs): how the 
phase manifests itself in an interference experiment (such as the one performed recently at 
MIT J5|), how the phase can be established by measurement, and how it evolves in presence of 
the elastic atomic interactions (see e.g. || and references therein). In this paper, in view of a 
possible experimental investigation of these problems, we complete the theoretical work already 
done on this subject by studying the dynamics of the relative phase in presence of loss processes 
occurring in the two condensates. Such loss processes, unavoidable in a real experiment, are 
due for example to collisions of condensed atoms with the background gas, or to three-body 
collisions between condensed atoms followed by recombination of two atoms to form a molecule 
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We consider two mutually non interacting and spatially non overlapping BECs in two trap- 
ping potentials. We suppose that the experimentalist has at hand a device, such as the one 
depicted in fig.[l], allowing both the measurement of the relative phase between the condensates 
and the preparation of a state with a well defined relative phase ||. Starting from an initial 
state with a well defined relative phase, we imagine that the two condensates evolve indepen- 
dently, under the influence of the atomic interactions, during a given time interval t at the end 
of which a measurement of the relative phase is performed. By repeating this procedure many 
times, one accesses the probability distribution of the relative phase ||. 

In the lossless case, the relative phase shows collapses and revivals periodically in time due 
to the presence of elastic atomic interactions. In presence of losses, we find that a few inelastic 
processes are sufficient to dramatically damp the revivals of the phase. In practice, for typical 
experimental configurations, the observability of the revivals is limited to condensates with a 
small number of atoms, of the order of a few hundreds, for which the revival time is of the order 
of 0.1 to 1 second. 



In section || we present the theoretical model describing the evolution of the system in 
presence of losses. An interesting feature of the model is that it can be solved almost exactly 
analytically within the Monte Carlo wave function approach recently put forward by several 
authors [|TI], [II], [T2|, |TB). We take advantage of this circumstance in the following sections, to 
deduce analytical expressions for the interesting phase-dependent measurable quantities, and 
to a give a simple picture of the phase dynamics in presence of losses: 

In section Q we find an approximate analytical expression for the evolution of a single 
stochastic wave function, and we give a simple physical interpretation of the result pointing out 
separately the role of the elastic atomic interactions and of the losses in the dynamics of the 
relative phase of the condensates. In sections f| and |5| we concentrate on the case in which the 
two condensates are placed in two identical traps and have initially the same average number 
of atoms, and we use the analytical results of section |3| to calculate the time dependence of 
some relative phase dependent quantities. In particular in section |] we consider an interference 
experiment where one counts the atoms detected in the two output channels of the beam-splitter 
of fig.|l|, and we analyze the two different physical situations in which the condensates' relative 
phase is initially sharply defined or is described by a "broad" relative phase distribution with a 
width ^> In section [| we imagine instead an experiment in which the time evolution of the 
whole relative phase probability distribution is measured. Sections |6] and [7] are dedicated to the 
analysis of additional features that would appear in an experiment; the effect of asymmetries in 
the parameters of the two condensates and in the initial average number of atoms is considered 
in section |6|, and the effect of fluctuations in the initial total number of atoms is considered in 
section [F[ Some concluding remarks are presented in section |S]. 



2 The model 
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2.1 The master equation 

Let us consider two mutually non-interacting and spatially non-overlapping BECs A and B in 
two harmonic potentials. Our starting point to describe the evolution of this system in presence 
of m-body losses is a master equation for the density matrix p describing the atoms in the traps: 



(i) 



where {X, Y} denotes the anticommutator, and [-^(f)]™ is the field operator raised to the power 
m which suppresses m particles in f. In second quantized form the Hamiltonian H reads: 



H = d 3 r 



(2) 



where Hq is the one-particle Hamiltonian including the trapping potential and the kinetic 
energy, and g = Anh 2 a/M where M is the mass of the atoms and a is the s-wave scattering 
length. 

The loss terms in Eq. ([!]) are parameterized by the number m of particles lost per collisional 
event and by the collisional constant k. Physically the case m = 1 corresponds to collisions of 
atoms in the condensate with atoms of background gas in the cell; the case m = 2 corresponds to 
spin-flip collisions between condensed atoms in magnetic traps, as only specific spin components 
are trapped; the case m = 3 corresponds to three-body collisions between condensed atoms, 
leading to the formation of an excited molecule and a hot atom supposed to leave the condensate. 
The collisional constant k for the processes m = 1 and m = 3 has been measured for 87 Rb atoms 
at JILA 0] and for 23 Na atoms at MIT ||. The collisional constant for the m = 2 process 
has not been accurately measured for these atoms yet, as the two-body losses seem to give a 
smaller contribution to the total decay rate. 

We assume that at any time the state of the condensate A (resp. B) can be described in 
terms of a single occupied mode, neglecting the excitations out of this mode due to a non-zero 
temperature or to the loss processes. We assume furthermore that these modes are the single 
particle ground state wave functions <p a , 0& given self-consistently as functions of the number of 
particles by the Gross-Pitaevskii equation: 

H + gN e \<f> e (r; N t )\ 2 } e (f; N e ) = p e {N e )<j> e {r- N e ) , (3) 

where the p e (N e )'s are the chemical potentials for the condensates with N e particles, and where 
the wave functions <fr € are normalized to unity. In more mathematical words we approximate 
the atomic field operator by: 

^(f)= E c e e (r;A>) (4) 

e=a,b 

where the operators c\ (cjj) and c a (q,) create and annihilate a particle in the condensate A 
(B) respectively, and where N e = c\c e are the operators giving the number of particles in each 
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condensate. Note that we keep in Eq.(f|) the dependence of the mode on the number of particles 
in the condensate. 

By substituting Eq.(|) into Eq.(Q) we get 

H = E a {N a ) + E b (N b ) (5) 

with 



E € (N € ) = N € 



d 3 r<j>* e (r; NjHoUr; N e ) + iV e )| 4 



(6) 



(we have used N e — 1 ~ N e ). 

By assuming that in the considered time interval the atom number distributions in the two 
condensates remain peaked around the initial average values: 

N e = Tr[p(0)ct Ce ] , (7) 

we expand the condensates' Hamiltonian around N a , N b keeping up to the quadratic terms: 

H(N a , N b ) ~ H q (N a , N b ) = J2 E (Ne) + (N e - N e ) fi e (N e ) + \(N e - N e ) V e (iV e ) . (8) 

e=a,b 

In our model we will use this quadratic version of the Hamiltonian, where the chemical potentials 
fi a and fi b and their derivatives can be calculated by solving numerically the Gross-Pitaevskii 
equation (||). 

We now substitute our ansatz Eq.(f|) in the loss part of the master equation; since the 
condensates do not overlap this amounts to the substitution 

imr - e ^ N e )r (9) 

e=a,b 

in Eq.(|]). In contrast to the Hamiltonian part which required a careful quadratization in 
N e — N e to get the correct phase dynamics, the dissipative part will be treated to lowest order 
by replacing N e by N e in Eq.(|9]). This allows us finally to obtain a master equation of the form: 

e=a,b 

where (for e = a,b) we have introduced the rates for the m-body collisions: 

lt = n ( d?r\U?\N e )\ 2m . (11) 
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2.2 Stochastic formulation 



To study the evolution of the system we adopt the Monte Carlo wave function point of view 
[ |IU| which provides us with a stochastic formulation of the master equation ( [X0|) . To this aim 
we introduce the jump operators: 

& = VtW e = a,b (12) 

and an effective Hamiltonian: 

H-^-YE S\S e . (13) 

e=a,b 

The Monte Carlo wave function \ip(t)) undergoes a non hermitian Hamiltonian evolution ruled 
by if e ff (plus a continuous renormalization) interrupted by random quantum jumps occurring 
at a rate (ip{t) \ J2 e =a,b(SlS e )\ijj(t)) , where \ip{t)) is normalized to unity The effect of a quantum 
jump is to replace by S t \ip) up to a normalization factor. Physically this corresponds to the 
loss of m particles in the condensate e via the m-body collisional processes described above. 
The two kinds of jumps e = a,b occur with relative probabilities: 

Pa = (m\sjs a \m) (u) 
p b m)\sis b m)) ' 

Starting with a state with a fixed total number of particles N, we can expand at each time the 
state vector on the Fock basis 

\m)= E d Na \N a ,N-N a ) , (15) 

N a =0,N 

where N is the total number of atoms at time t in the two condensates, and we can carry 
out the evolution numerically. The mean value of an observable O is obtained by averaging 
the expectation value (ip(t)\0\tp(t)) over all possible stochastic realizations for the evolution of 

Usually the Monte Carlo wave function technique is carried out purely numerically. It turns 
out that for the present problem it is possible to treat analytically the evolution of a Monte 
Carlo wave function and, after a minor approximation, average analytically over all the possible 
stochastic realizations. This leads to a simple interpretation of the dynamics and allows the 
derivation of analytical formulas for observables' mean values. As it will appear in the figures 
the analytical results are in good agreement with the numerical results. 

3 Evolution of a single wave function 

In this section we derive an approximate formula for the evolution of a single stochastic wave 
function, and we discuss its physical interpretation. We first consider the simple case in which 
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the condensates are initially in a phase state, introduced in the beginning of the section, and 
subsequently the general case in which the initial state is characterized by a given relative phase 
distribution. 

For the following it will be useful to introduce the operators 

N = N b + N a n = N b -N a (16) 

corresponding to the sum and difference of the number of atoms in A and in B. 



3.1 Phase states 

A very useful class of states of two condensates is represented by the phase states [ |14| : 



l^=7Pff (4^ + cIe-*T|0> (17) 

having a fixed total number of particles N and leading to a well defined relative phase 2</> 
between the condensates A and B. These states have the remarkable properties: 

c^)n = ^e^^- s ^\^) N ^ e = a,b (18) 

e~ iah \4>) N = \4> + a) N Vor, (19) 

where the 5 e ^ for e, e' = a,b are Kronecker deltas. The first property reflects the fact that in a 



phase state, all the particles are in the same state (see Eq. (|17D ), and the second one shows that 
n and <fi are to some extent conjugate variables like the momentum and position of a particle. 
Note that the phase states are not orthogonal: 

N ( ( p'\4>) N = [cos(<f>-<f>')) N , (20) 

though the function [cos(0 — <f)')] in Eq.(|20|) becomes very peaked around zero when N — ► oo 
with a width scaling as 1/y/N. Any state with a total number N of particles can be expanded 
on the overcomplete set of phase states: 

l-Ti/ 

where c(4>) can be obtained from the expansion of the state vector on the Fock state basis: 

c^)=A- 1 £ 2 AV2 f ^N -N a )\ \ 1/2 et{N _ 2Na) , {NafN _ NM (22) 



a r /2 ^ C (d>) \<i>) N , (21) 

J-tt/2 7T 



N a =0,N 



The quantity |c(0)| 2 can be interpreted as the relative phase probability distribution 0. This 
distribution, flat for a Fock state and very peaked for a phase state, is normalized in such a 
way that: 

J" 2 , * M#)|" = 1 . (23) 



The factor A in Eq.(|21|) ensures that \ip) is normalized to unity. For iV 3> 1 and for a c(0) 
varying slowly at the scale l/\fN, we can replace the scalar product n{4>'\4 > )n by the delta 

distribution ^J^m/N8{(j) - 0') to obtain A = (^f^- 
3.2 Approximate expression for \4>(t)) 

Consider the evolution of the state vector \tp(t)), from a time to = to a time t, for a particular 
stochastic realization. We imagine that k quantum jumps, each corresponding to the loss of m 
particles, occur at times ti, ...,tk separated by time intervals Tj = tj — with j = 1, k; the 
k th jump takes place in the condensate e& with e& = a, b. We have: 

= Me-fc^-^S^e-te^S^e-* 11 ^-!... S ei e-^ H ^\^(0)) (24) 

where A/" is a normalization factor. By using the identity: 

[c e ] m f(N a , N b ) = f(N a + m5 £ , a , N b + m8 e , a ) [c e } m e = a, b , (25) 

we shift all the jump operators in Eq.(|24D to the right by letting them "pass through" the 
exponentials and we obtain: 

= Afexp[-iH eS ({N e })(t - t h )/h) ex P [-iH efi ({N e + mS e , ek })T k /h) 

exp[-iH eS ({N e + m(8 e , ek + S^^n-i/h] ... ]J S ej \^(0)) . (26) 

j=i,k 

We introduce now the major approximation in our calculations by replacing [cj] m [c e ] m by N™ 
in the expression for the effective Hamiltonian Eq. flT5|) , supposing that the fraction of lost 
particles is small. The resulting effective Hamiltonian then takes the form: 

H eS = H q - y A , (27) 

quadratic in N a and Nb, where A is a constant representing the mean total number of collisional 
events per unit of time: 

A = A a + A b with A a = 7a iV a m , \ b = lb N b m . (28) 

In this approximation the statistics of the quantum jumps is simply Poissonian with a parameter 
A and 5 b>e , = 1 — 6 a>ej takes the values 1 and with probabilities X b /X and A a /A respectively, 
according to Eq. ([14]) . 

We then expand the effective Hamiltonians in each exponential in Eq.(|2"6D around A^ a , N b 
in powers of m5 etek , m(5 e ^ k + 5 £)efe l ), etc. Due to the quadratic dependence of Eq.(p7[) on N a 
and N b we limit the expansion at the first order, the subsequent terms being constants or zero. 
By using Eq.(^) we then obtain the following result for the state vector at time t: 

m))=xe- xt/2 m)Ui{t) n ^#(°)> • (29) 

j=i,k 
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In Eq.fl29|) we have introduced the unitary operators 

U (t) = exp[-iH«({N e })t/h] (30) 

(31) 



Ui(t) = exp 
where for e = a, b: 



A e = m 6 ^l T 3 = m J2 5 e,eiU (32) 

j=l,kl=j,k l=l,k 

are random quantities that depend on the particular realization. 

We sketch out briefly the physical interpretation of the result Eg. (|29"D , considering the action 
of the successive factors in Eq.(p9|) on a phase state defined in Eq.([l7|). 



The factor Uo(t) in Eq. (j29|) accounts for the evolution in absence of losses. Expressed in 
terms of the operators N and n of Eq. ([TBI) it involves: 



H q ({N e }) = f (N) + hv{N) + n\^ h + //J/8 . (33) 
We have used Eq. (|]) and we have defined 

v(N) = ^{fi b ^^(N -N)- ^±^(N b - N a )} , (34) 

where N = N a + Nf, and where /x e stands for fi t (N t ). From the properties of the phase 
state we find that the terms in n and n 2 in Eq . (j33|) , when exponentiated in Uq, (i) shift 
the relative phase at the iV-dependent constant speed v(N) and (ii) spread the relative 
phase (in a way analogous to the spreading of a wave packet of a massive particle under 
free evolution), respectively. The term fo(N) in Eq. (|33]) is a function of the total number 
of atoms iV only and plays no role, since it amounts in Uo{t) to adding a global phase 
factor to the wave function. The phase- spreading will eventually lead to a collapse of 
the relative phase [Q. On the other hand due to the discreteness of the spectrum of 
the operator n (the spectrum of n consists of even integers for an even N, and of odd 
integers for an odd N), there are special times at which the exponential operator Eq.(|3"3"D 
reduces to a mere translation of the relative phase, yielding the well known result that 
revivals should follow the collapses of the relative phase. More precisely if one uses the 
expansion Eq.(|T5|) for the phase state defined in Eq. (|T7|) , one realizes that a relative phase 
distribution initially peaked around O displays revivals at the times: 

tR = qn/x, 1 integer (35) 

where we have introduced: f / 

X = —2h-- (36) 



8 



At these times, for N even: 



and for iV odd: 



-ixn 2 t R /A\ 



-ixn 2 t R /A\ 



) N = |0 + g7r/2) 



N 



(37) 



e — ' 1WJV = e-^ 4 |0)iv (38) 

The initial relative phase distribution is then reconstructed around (0o +v(N)tR + qn/2) 
for N even and around (0 O + v(N)tn) for iV odd. 

The factor £/i(i) in Eq.(^) accounts for the presence of losses. Expressed in terms of the 
operators n and N, it involves: 



dH q dH q 

-({N e }) A a /h + jr-m}) A b /h = J\(N) + nD 



dN a 



dN h 



(39) 



where global phase factors are included in f\(N). The translation operator n appears in 
Eq.(E^) multiplied by a random quantity D defined as: 



D = m ti 

l=X,k 



2h 



(40) 



Equations ( p!9|) and (|39|) show that the relative phase in a single stochastic realization is 
shifted by the random amount D due to the loss processes. This effect will turn out to 
have a dramatic influence on the coherence properties of the condensates. 



Finally in Eq.(p9|) the action of the jump operators on a phase state is simply: 



n s, 

j=i,k 



NN-l N-mk + l 



2 2 

where we have introduced the quantity 



1/2 



>N-mk 



a = m [ 2 $b,e 
j=i,k 



1 



(41) 



(42) 



Apart from numerical factors that will be absorbed in the normalization and the phase 
factor involving a, Eq.(|4l|) amounts to reducing by a random amount the total number 
of particles. 

In the general case, an initial state with iV particles can be expanded on the phase states set 
(see Eq.([H])). By using Eqs. (|33|) , (|39|) , (|4i"D , and getting rid of the global phase factors we 
then obtain the wave function: 



| m) = B(t) f /2 ^c(0, 0)e-^ 2 */ 4 e-^|0 + D + v(N - mk)t) N - mk 

J-tt/2 7T 



(43) 



where B(t) is a normalization factor. 
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4 Mean beating intensity of the condensates 



To monitor the evolution of the relative phase between the condensates, a possible choice is to 
determine the relative phase dependent quantity (<4c&) after some time during which the two 
condensates, initially prepared in a state with a defined relative phase, evolve independently. 
As the relative phase between the condensates is affected by the elastic atomic interactions, the 
average (c^q,) undergoes collapses and revivals in time. 

In the situation described in fig.|l| the measure of (c^q,) would correspond to the following 
measurement scheme: Prepare a state in which A and B have a well defined relative phase 
[f|; let the condensates evolve during a time interval t; then let p -C N atoms escape from 
the condensates and beat them on the beam-splitter. The counts registered in the two output 
channels of the beam-splitter will be fluctuating variables whose averages over many realizations 
of the whole procedure are 0: 

,V (4±4)( C a±Q>) pi 



jV 2 N2 

The difference between J + and J_ gives then the real part of (c\cb). 



) * "bo + (cfa) ± 2Re(ct Cb ) , (44) 



We shall now use the approximated formulas ( p5|) and to calculate the time dependence 
of (<4q,). The main result of this section is that the revivals in this quantity are damped in 
time with a simple exponential law e~ xt where the constant A, defined in Eq.(^), is the mean 
number of loss processes per unit of time. 

In the present and in the following section we restrict for simplicity to the perfectly symmet- 
ric case where the two trapping potentials are identical and the two condensates have initially 
the same mean number of particles: 

N a = N b , (45) 

la = Jb, (46) 

Ha = /J, b . (47) 

Moreover we consider an initial state having a fixed total number of particles equal to iV; and 
as a reminder of this choice (when it is the case) we will attach a superscript (...) fix to the 
averages. The non symmetric case for the condensates will be considered in section ^ while 
the effect of fluctuations in the initial total number of atoms (requiring a further averaging over 
N) will be analyzed in section [7|. 

We calculate (cjjq,) fix in two different physical situations. The first one refers to a sharply 
defined initial relative phase (A<p — for which we choose a phase state as the initial state; 
the second one, probably more realistic from the experimental point of view, makes use of an 
initial phase distribution much broader than In each case we first calculate the expectation 

value of the operator O = c\c\, for a single stochastic realization using the results of section 
and then take the average over the stochastic realizations. In the whole paper we will 
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denote with (i[)(t)\0\i/;(t)) the single realization expectation value and with (O) the quantum 
mechanical average. 



4.1 Case of an initial phase state 

Let us assume \i/j(0)) = \</>)n\ by using equations (|29"D and (|33D, (j39|), (f4l|), for a single realiza- 
tion, we find: 

(m\c\c b \m) = N- m k^ + D\e^ h2t cic b e-^ flH \ ( l ) + D) N ^ mk (48) 

where \ and D are defined in Eq.(^) and Eq.(^) respectively. Note that the contribution 
involving the drift velocity of Eq. (|3^) vanishes as we are considering here the symmetric case. 
The quadratic dependence on n in Eq.(^) can be eliminated by shifting c\c b through the 
exponential e - **' 1 1 using Eq.([25|): 

e^clcbe-^ = e- ix{fl+1)t clc b (49) 

so that 

(m\c j a c b \m) = N-mk((p + D\e- ix{h+1)t c j a c b \(p + D) N _ mk ; (50) 
by using the properties (|18|), (|T9|) and (^Cj) we then have: 

m)\clc b m)) = e~^e- 2W [ C os( X t)] N - mk - 1 . (51) 



The next step is to take the average of the result Eq. (|5T|) over the stochastic realizations which 
amounts to averaging over the random variables k, tj and 8 b ^ j (the last two variables appearing 
in the random quantity D). We show the calculation of the average in detail in the appendix 
A. The result for (<4o,) fix reads: 

{c\c b r = e-^e~* £ «(*)]* [cos^)]^ 1 , (52) 

k=0,N/m-l / K - 



where the function u(t) is given by: 



u(t) = . (53) 



By identifying the factor iV — mk with N under the assumption of a small fraction of lost 
particles, and by extending the sum over fc up to oo, we are able to express the result in a 
compact way: [] 

( c t Cfo )fix = e -2i<P e -Xt[l-u(t)/cos^ x t)]E[ C0S ^ x t)f-l ( 54 ) 



1 It should be noted however that the compact formula (|5J) diverges for \t = n/2 + qn, where the explicit 
sum Eq.(52) should be used instead. At such points (c^Cb)^* = anyway. 
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The factor [cos(xt)] 7V_1 in Eq.(|54]), already obtained in [15| in the absence of losses, is respon- 
sible for the collapses of the average value {c[cb) &x and for revivals at times = qn/x with q 
integer. The collapses and revivals of (c^q,) are shown in fig.0 both (a) in absence and (b) 
in presence of three-body losses. We see immediately that the losses have a dramatic effect 
reducing exponentially in time the average with the rate A given by Eq. (p8|) . In fact at a re- 
vival times t = tn, u(t) vanishes so that the average value of (ctq,) fix is simply attenuated with 
respect to the lossless case: 

(c\c b )t tR = (-l) q{N ~ 1] (clc b )t e~ xt « , (55) 

by an exponential factor which is exactly the probability that no particles are lost up to time 
t. The dramatic effect of losses on the revivals, already when Xtn ~ 1 (that is one loss process 
has occurred on average at the revival time), can be understood by the fact that in each single 
realization experiencing a quantum jump at a time t ~ the relative phase is shifted by an 
amount D ~ 7T. This point will be further exemplified in section H. 



4.2 Case of an initial relative phase distribution broader than that 
of a phase state 

Since it may be difficult to prepare experimentally the condensates in a phase state we now 
consider the more realistic case in which the initial relative phase distribution |c(0, 0)| 2 for the 
condensates is broad as compared to 1/yN. To be specific we assume that the initial relative 
phase distribution is a Gaussian centered at <fi = 0: 



c(0,O) = £ exp (-0 2 /(4A0 2 ) 



< Acf) < 1 



(56) 



where ranges between —ir/2 and tt/2. This choice corresponds to a Gaussian distribution for 
the number of particles in the condensates: 



(iV a ,iV-iV a |^(0)) = ^e-( Ar - 2 ^) 2 / 4A " 2 



(57) 



with An = 1/2. 

For a single realization, we use Eq.(^3|) and we proceed along the lines of the previous 
calculation to get: 



(m\c{c b \m) 



~irN~ 


1/2 


2 





T/2 WW c y Q y {M N-MW 

tt/2 71 71 2 



(5f 



where N = N — mk with k equal to the number of quantum jumps experienced by the Monte 
Carlo wave function up to time t. Now by using the fact that the scalar product between the 
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phase states for N 3> 1 is a very peaked function of — 0' with respect to the other functions 
in the integral, we perform the substitution: 



2n 

n-i{4>' - xA4>)n-i -> cos^^govr) \j -^${4> + Qo^ ~ X* ~ 0) ( 59 ) 



where the integer q is chosen such that — n/2 < (%t + — q n) < n/2. As the factor c(0, 0) 
defined in Eq . (|56|) is peaked around = 0, we neglect the dependence of q on so that the 
integer q is finally chosen such that —tt/2 < (xt — go 71 ") < 7r/2. In this way we obtain 

m)\c\c b m)) = (-i)^-i)^ e *(x«t-2D) 

r' 2 !tt c U o) c * (0 + x t - go7r , 0) e-^xt-go*) . (60) 

J-tt/2 7T 

The next step is to average the factor e^ xat ~ 2D ^ over the stochastic realizations. The procedure 
closely follows the one in the appendix A. By identifying N with N, as in the previous case, 
and by extending the boundaries of integration in Eq. (|60D to ±oo we can express the result in 
the compact formF]: 

AT +°° 

( c t Cfe ) fix = —e^ 1 -"^)] ST e -[(x*-9T)/2] 2 /2A^^_ 1 ^(iV-l) ^) 
" 2 9=0 

where u(t) is defined in Eq. fl53|) . The factor involving the sum over q in Eq.(|62|) plays the role 
of the factor [cos(xt)] N ~ l in Eq.(|54]) which was obtained for an initial phase state. At each time 
tfi = qn/x there is a revival of the quantity (c] l q ) ) fix and Eq.flB^D reduces to the very simple 
expression: 

(c[c b )t tR = (-ly^Hclc^to e~ xt « . (63) 

This formula does not depend on the initial width A0 and coincides with the one Eq.Q55D 
obtained for a phase state. There is therefore no possible way of reducing the damping of the 
revivals by adjusting the initial width of the phase distribution. Only the temporal width of 
the revivals is larger for a distribution broader than that for a phase state, as it clearly appears 
from a comparison between fig.|] and the previous fig.@b. 



Remark: Formula (|62|) can also be used to study the collapse of the phase around t — 0. 
For short times (t <C t#) we expand u(t) to second order in t obtaining: 



2 1 1 8A0 2 



3 



(64) 



To obtain Eq.(p2) we use the condition A(j) << 1 to set: 



(«I&4(£;^„,0) e -)4. (6!) 
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In the absence of losses we recover the collapse time t c = 2A<p/x @ as the half temporal width 
at the relative height e" 1 / 2 of the mean signal (c^cj,) flx . Losses start accelerating the collapse 
significantly when m 2 A(f) 2 Xt c > 1. In this regime of course the subsequent revivals cannot be 
observed. 



5 Evolution of the relative phase distribution 

We turn now our attention to the phase distribution |c(0)| 2 which could be reconstructed in 
an experiment for example via a series of multichannel measurements. We show an example of 
the procedure in fig.|] |16|], 0. 

In the frame of our model, the evolution of c(<ft) can be obtained numerically from the 



evolution of the state vector \ip(t)) expanded on the Fock state basis by using Eq. (|22|) ; however, 
as we show in the following, the approximated analytical treatment allows us also in this case 
to find some simple results at the revival times. 

Let the initial state of the condensate, with a total number N of atoms, be characterized 
by a given relative phase distribution c(0, 0); the state vector at time t is then given by our 
approximated formula Eq. (f£5|) . One can easily check that the integrand in Eq.(|4"5|) is periodic 
of period 7r so that we can shift the interval of integration to obtain: | 

\i,{t)) = B{t)e-^ 2t " f 72 — c(0 — D — v(N)t, 0)\<f>) s (65) 

J-w/2 7T 

where c(0) = e~ ia ^c((j)) and N = N — mk. This result has a very suggestive interpretation: 
The loss processes in a single realization shift the relative phase distribution by a random 
amount D, and the overall evolution can be separated in a random shift plus the Hamiltonian 
evolution. To make clearer this interpretation, we have plotted in fig.[5] the phase distribution at 



the second revival time (given by Eq. (|35] ) with q = 2) for different realizations. For Xtfi ~ 1, as 
in the figure, there is an important fraction of realizations in which the relative phase is shifted 
considerably. This is the reason why the relative phase distribution at the revival time will be 
smeared out by the losses when we take the average over the stochastic realization, which we 
do now. 

As in section f| we consider the symmetric case defined by the Eqs.([EJ), (f46"D, (f4?D. Fur- 
thermore we restrict ourselves to the revival times t = tn = qn/x, Q integer (see Eq. (|35|)). 
In this case the Hamiltonian evolution operator in Eq.(^) takes a simple numerical form (see 
Eq.(|37|) and Eq. (p8|)) and by comparing Eq. (|65"D to Eq.(^l|) we can simply read out the phase 
distribution amplitude c((j),t): 

c{<t>,t R ) = 5(0^ -AO), (66) 



When — > + 7r, c(4>, 0) is multiplied by (— 1) , exp(— ia4>) is multiplied by (— l) mfe , and the phase state 

JV 



\4> + D + vt)ft is multiplied by (-i)N-mk 
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where: 



0jv — — ?7r/2 for A" even (67) 
0^ = for N odd. (68) 

From Eq.(pB|) we see again that a single loss event (which can lead to D ~ 7r) has a dramatic 
effect on the phase distribution. 

As shown in the appendix B the phase distribution at the revival times averaged over the 
stochastic realizations takes the very simple form: 

(\c^,t R )\ 2 f x = (l-e xt «) + e- XtR \c^ N ,0)\ 2 . (69) 

At the revival time the relative phase distribution is "damped" by the factor e~ xtR while a flat 
background component appears. This effect is clearly shown in fig|], where we have plot the 
averaged relative phase distribution at t — and at the second revival time. 



6 Effect of an asymmetry between the two condensates 

In the previous sections we have investigated the relative phase dynamics in the symmetric case 
for the two condensates. In this section we extend the analysis to account for a small imbalance 
in the initial average number of particles 

\N b -N a \ <AT, (70) 

where N is the average of the total initial number of particles, and for arbitrary values of the 
parameters /i a , fib, 7a, 7b- We restrict the calculation to the contrast of the interference fringes 
between the two condensates averaged over many experimental realizations, assuming an initial 
phase distribution broader than the phase state. 

Our initial Monte Carlo wave function has a fixed total number of particles equal to N, and 
a Gaussian distribution for number of particles in each condensate. The calculation of (c^q,) 
is now slightly more involved than in the symmetric CcLSG, clS the phase distribution amplitude 
c(0, 0) acquires a phase factor varying rapidly with at the scale 1/y/N. All the calculations 
are therefore put in the appendix C, and we give here the result only at the revival time t = t R : 

{c\c b )tt R = (-l)W-V-e-™Wn e -xt a [i-uit a )] ; (71) 

where v(N) is defined by Eq.(|34l) and U{t) is a function of time (see Eq.(|9~7D in appendix C). In 
fig.0 we show an example of the time evolution of (cjjC ) fix in the case of a 10% asymmetry in 
the initial number of particles N a and N^. As far as the damping of the revivals is concerned, no 
significant difference appears with respect to the symmetric case. The damping of the revivals 
is in this case ruled by the exponent: 

- Xt R [l - ReU(t R )} (72) 
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where: ^ 

ReU(t R ) = - ( A& smc(mfi' b t R /h) + X a sine (m//^/ ft) ) , (73) 

where sinc(a;) = sin(x)/x. Obviously |Ref/(t^)| < 1, meaning that an asymmetry between the 
condensates cannot amplify the revivals with respect to the lossless case. From Eq. ( |73|) we 
notice, just as a curiosity, that a complete suppression of the effect of the losses (ReU(t R ) = 1) 
would occur only in the case in which there are no losses in the condensate A (A a = 0) and no 
elastic interactions in the condensate B (fi' b = 0) (or vice versa). 

A trivial effect of the asymmetry, evident in fig.[?|, is the appearance of oscillations of the 
mean value (c^Cfe) fix due to the non zero drift velocity of the relative phase of the condensates. 
We will see in the next section that this effect, harmless at first sight, can have dramatic 
consequences when we consider the effect of the dispersion in the initial total number of particles 
N. 



7 Effect of fluctuations in the total number of particles 

Through all the previous sections in this paper we have chosen an initial state, represented by 
our initial Monte Carlo wave function, with a fixed total number of particles in the condensates. 
The averages that we calculated (...) fix then correspond to the real quantum mechanical averages 
supposing that the initial total number of atoms is fixed to a value N for any realization of the 
experiment. In practice it is probably difficult to produce a Fock state for the condensates and 
the total number of atoms should be governed by some probability distribution P(N). Since 
we have analytical formulas for the quantities of interest (such as the average (c^q,)^), it is 
very simple to add a further averaging over iV for a given P(N). Suppose for example that 
the distribution for the initial total number of atoms is a Poissonian distribution of parameter 
N. By averaging the result Eq.([TI]), valid at the revival times t R for slightly asymmetric 
condensates, we get: 

|( c t c ^Pojss| = ^_ e -\t a [l~RcU(t R )] e -7V{sin 2 (^t ii /2ft)+sm 2 (^t fl /2a)} _ ^ 

The result Eq.(^) shows that a slight asymmetry between the condensates kills the revivals 
of (c^Cft). This is due to the fact that the drift velocity of the relative phase v(N) in Eq.(|7l|) 
depends on the initial total number of particles, giving to (c^q,)^ a phase factor of the form: 

exp[-2iv(N)t R ] oc exp[* (N - N) t R ] = exp[* (N - N) 4^4 ^ ■ ( 75 ) 

To be able to observe the revivals it is then necessary to be as close as possible to the symmetric 
conditions in order to satisfy: 

4— < 1 , (76) 
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where AN is the width of the distribution P(N). 

If the symmetry between the condensates is perfectly realized, the atom number fluctuations 
have the simple effect of doubling the revival time. We show an example in figj^ where we 
averaged the result for (c^q,) 6 * for an initial phase state (fig.||) using a Poissonian distribution 
for P(N). The main effect is the disappearance of the "odd" revivals; this is due to the fact that 
the amplitude of these odd revivals for N particles is proportional to [cos(g7r)]^ Ar_1 ' ) = (— l)^ -1 ) 
which alternates its sign depending on the parity of N. 

In fact it is possible to show that a Poissonian ensemble of phase states is equivalent to a 
coherent state for the two condensates, as long as one calculates the mean values of operators 
commuting with the total number of particles in the condensates. For the perfectly symmetric 
case in fig.|]we then recover the result obtained in |1| (in the absence of losses) i.e. the doubling 
of the revival period for a coherent state of the condensates as compared to the phase state. 



Within the coherent states pictures we can also reinterpret the result Eq.(74) for the asym- 
metric case in the following way: in order to observe a revival of the relative phase between 
two condensates it is necessary that both condensates display a phase revival at the same time 
i.e. (j,' a /2htn = qn and fi' b /2htR = q'n, with q, q' integers. 



8 Concluding remarks 

We have studied the dynamics of the relative phase between two Bose-Einstein condensates in 
presence of m-body loss processes in order to question the observability of the collapses and 
revivals of the phase predicted by purely Hamiltonian models. 

We have shown that the losses damp exponentially in time the phase dependent quantity 
(c\cb) (see Eq.(j55|) for an initial phase state and Eq . (|63|) for an initially broader phase distribu- 
tion). The decay rate A of (c^q,) coincides (up to the factor m) with the mean total number of 
particles lost per unit of time, and it is therefore approximately TV times larger than the inverse 
lifetime of a particle in the condensates, where N is the total number of particles initially in 
the condensates. 

The dramatic effect of the losses on the relative phase has been suggestively interpreted 
within the Monte Carlo wave function approach. In a single realization each single loss event 
occurring at a time of the order of the revival time shifts the relative phase by a random 
amount of the order of n. A few loss processes are then sufficient to smear out the relative 
phase completely at the revival time when the average over the stochastic realizations is taken. 
For this reason the experimental observation of the revivals is limited to condensates with a 
small number of atoms where the condition < 1 (where tu is the revival time Eq.(|35|)) can 
be satisfied for all the relevant loss processes in the system. 

In order to give an idea of the possible scenarios and of the order of magnitudes in different 
experimental conditions, we have shown in fig.^| the loss rates due to one-body and three- 
body collisions and the inverse revival time as functions of the total number of atoms, for 
two different values of the trap frequencies. For higher trap frequencies (fig.|9]a) the revivals 
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occur on a shorter time scale and one is confronted mainly to three-body losses, while for 
less confining traps (fig.|9|b) collisions with the residual gas should be taken into account due 
to longer revival times. Fig.^J] shows that phase revivals in presence of losses are in principle 
observable in condensates with some hundreds of atoms. 

By studying the general case of two asymmetric condensates, and the effects of fluctuations 
in the initial total number of atoms in the condensates, we have finally pointed out a practical 
difficulty which should be overcome in order to observe phase revivals. The difficulty comes 
from the fact that in the case of two non perfectly symmetric condensates their relative phase 
drifts with a velocity depending on the initial total number of atoms. For this reason random 
fluctuations in the initial number of atoms turn out to destroy the relative phase revivals when 
the asymmetry is too large. A possible way to overcome this problem is of course to use two 
almost symmetric condensates. Another possibility, which we have not examined in detail, 
would be to use a condensate A which has a collapse time longer than the duration of the 
experiment (N(fi' a tn/h) 2 <C 1) as a phase reference to measure the evolving phase of the other 
condensate B. 
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Appendix A: average of the phase factor e 



-2iD 



In this appendix we derive the average over the stochastic realizations of the quantity e~ 2tD S(k) 
where D is defined in Eq. fl40|) and where S(k) is an arbitrary function of the number of jumps 
k. We perform the average over the variables 5& i£ . first, using their probability distribution 
given after Eq . (^8|) ; we have: 

(e~ 2iD U,e, = II \{\ b e~* m ^ + \ a e* m ^) = ]J f(tj) . (77) 

3 j=l,k * j=l,k 

In order to perform the average over the variables k and Tj, we need the probability distribution 
Pt(k, ti,t2, ■■■tk) of having in the time interval (0, t) exactly k jumps separated by time intervals 
Tj = tj — tj-±. Since we assume that the loss processes occur randomly with a constant rate 
A, corresponding to a waiting-time distribution of the form w(r) = Ae _Ar , the probability 
distribution P t (k, t±, t 2 , ■■■tk) is simply fl"3|j : 



P t (k,t 1 ,t 2 ,....t k ) = \ k e- Xt . (78) 
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Using this result we are led to the calculation of a multiple integral of the form: 



o<t 1 <t 2 ...<t k <t 



f(t 1 )f(t 2 )...f{t k ) dt x dt 2 ...dt k 



(79) 



where f{t) is the argument of the product in Eq. (|77D . Since / is equal to I a calculated for any 
permutation ...to-(fc) of the integration variables, we can write it as: 



We then obtain 



E / f(ti)f(t 2 )...f(t k ) dhdt 2 ...dt k 

a Jo<t aW <...<t a(k) <t 



(S{k)e^ D ) k , T]A =Y.S{k)^\f f{t) dt 



1 

hi 



f(t) dt 



(80) 



-Af 



In this last equation we may have to introduce by hand a cut-off N/ m — 1 over the index k if 
S(k) has divergences for k > N/m (i.e. when no particles are left in the condensates). 



Appendix B: phase distribution at revival times 

We are interested in calculating the phase distribution at the revival time averaged over the re- 
alizations that is (\c(<fi, tji)\ 2 )k,Tj,s b e . • We restrict to the symmetric case between the condensates 
and we start from Eq. fl66|) . By using Eq. (p2|) for t = we have: 



(\c(<P,t R )\ 2 ) k , Tj , 5 , 



\A0)\- 2 E E fac(iV a )fac*(iV:)(e 2j ^-^)^- D )) / 

N a =0,N N'=0,N 



where we have introduced the notation 



1 /2 

fec(JV a ) = 2 N ' 2 ( NJiN ~ Na)l ) (N a , N - N a \m) 



(82) 



(83) 



The calculation of the average over the stochastic realizations closely resembles the previous 
one Eq . (|77|) that we have explained in the appendix A; we have: 

k 



)fc,r,- 



E< 

fc>0 



-Mr ( Xt R. 
k\ 



sin[(A^ - N a )m X t R ] 
{N' a - N a )m X t 



R 



54) 



We note that the terms in the sum in Eq.(R3) for k ^ are equal to zero unless (A^ — N a ) = 
in which case the average in Eq.(|84|) is equal to one. We can then rewrite the result (^) as: 



(\c(M\ 2 ) 



2\fbc 



1-4(0)1 



E E ^,7vJfac(iV a ) 

N a =0,N N^=0,N 



+ (1 - 6n L ,nJ (fac(iV a )[fac(iV:)]*e 2 ^-^ 



(85) 
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Now by using the property: 

E |fac(iV a )|V(0)|- 2 = l (86) 

N a =0,N 

coming from the normalization condition Eq.(|23|) and from Eq. (|22|) , we find the suggestive 
result Eq.(H). 



Appendix C: asymmetric condensates 

In this appendix we show the explicit calculation of the mean contrast of the interference 
fringes (cJjQ>) fix for asymmetric condensates. We consider an initial Monte Carlo wave function 
for which the total number of particles iV is fixed and the number of particles in condensate A 
has a Gaussian probability distribution: 

(N a , N - AgV(O)) = g e -^~ x ^ 2 / An2 (87) 

where Q is a normalization factor and An is the standard deviation for the difference n in 
the number of particles in the two condensates. The quantities x a = N a /(N a + N b ) and 
Xb = Nb/(N a + Nb) are the average fractions of particles initially in the condensate A and B 
respectively, which are supposed to be fixed from one realization to the other even in presence 
of fluctuations of the initial total number of atoms. 
We suppose in what follows that 

1 < An < VN , (88) 

and 

\x a N - x b N\ < N . (89) 



We first derive the phase distribution amplitude corresponding to the initial state Eq.(| 
by using Eq. ([22]) . We evaluate the factorials in Eq. ([22]) using the Stirling's formula, and we use 
a local approximation valid for \N a — x a N\ <C y/N : 

NJ{N-N a )\ ^ (x a N)\(x b N)\ (Na _ XaN) Hxa/Xb) 



By approximating the discrete sum in Eq. ([22]) with an integral over N a ranging from — oo to 
+oo we obtain: 

c (0, 0) = Me-^^e^ (91) 

where: ^ 

k= (x b - x a )N - -^An 2 ln(x a /x 6 ) (92) 

and where M is a normalization factor obtained from Eq. (|23|) . We note that in the symmetric 
case Aa = N b , we recover the Gaussian dependence for c(0) of Eq.(^) with AnA<p = 1/2. 
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We are now ready to calculate {c\cb) Kx starting from Eq. fl43|) . The calculation closely follows 
the one in section f|. In particular we use the key property Eq. (|49f) to obtain: 



7T 2 ' ' ' ' J-tt/2 J-tv/2 



7! 1 - d( Ll e -{<t> 2 +<t>' 2 )^ e i{K- a )(<t>-<t>') 



(m\cic b \m) = -m)\'w\- 

| e^^^^^-xt^h-v (93) 



The phase factor e lK ^~^'^ in the integrand varies rapidly with — 0' at the scale 1/V~N when 
Nb — N a is larger than \^N. For this reason we approximate the scalar product between the 
phase states \<j>) $ and \(j>') ^ by a Gaussian exp(— iV(0 — 0') 2 /2) rather than by the 5 distribution 
of section HI. This leads to the approximation 



N-lW ~ Xt\<P)ff-l ~ (-l)^-l) e ~(N-lW-^ X t +qo nf/2 (Q4) 



where the integer go is chosen such that —n/2 < (xt — go 71 ") < it/ 2. By extending the limits of 
integration over 0, 0' to ±oo in Eq. (|93| ) we are then left with a double Gaussian integral that 
can be calculated exactly. The result is quite involved but it can be simplified by using the 
condition fl89| ) and Eq. (|88|) . We take the average over the stochastic realizations and we use 
again Eq.(jS8[) to simplify the result. We calculate the normalization factor B(t): 

1 ~ - 2 \Af\ 2 \B\\t) f^LV /2 f_i^y /2 e-§(«-) a /(^ + An") . (95) 
7T 2 ' 1 1 1 w V4An 2 / \N + An 2 ) V ; 

We finally obtain for the mean contrast of the interference fringes between A and B as: 

+oo 

(4c 6 ) fix ~ e -\t e -2iv(N)t y2 e -±An 2 [(xt-w)] 2 (_iy{N-i) 

q=0 

ES" 1 ^e~^ xt ~ qn) -^^[XtU(t)} k (96) 
where the function U(t) is given by: 

I / e imn' b t/n _ j e -imp' a t/% _ ]\ 

m = a ( A '^^7T + x " - mM n ) ' (97) 
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Figures 



Figure 1: Two BECs A and B in two non overlapping trapping potentials. Some atoms can be 
let out of the condensates towards a 50-50 atomic beam splitter. The detection of the atoms in 
the output channels of the beam-splitter realizes a measurement of the relative phase between 
the condensates. 



Figure 2: Collapses and revivals of (c* a Cb) for an initial phase state (a) without losses and (b) 
in presence of 3-body losses. The calculation is performed for 87 Rb atoms in the F — 1, mp = —1 
state and for isotropic harmonic traps. The 3-body loss rate is inferred from the experimental 
data of JILA. The initial total number of atoms is N = 301, and the harmonic frequencies 
are Q a /2n = Q b /2n = 500Hz. Diamonds: numerical result with 2.5 x 10 4 Monte Carlo wave 
functions. Solid line: analytical result. 



Figure 3: Collapses and revivals of (c^Cb) fix for an initial phase distribution broader than that 
of the phase state. The initial total number of atoms is N = 301. The initial distribution for 
the difference in the number of particles in the two condensates is Gaussian with a standard 
deviation An = 6 and a vanishing mean (so that N a = Nj,). The other parameters are as in 
fig^b. Diamonds: numerical result with 2.5 x 10 4 Monte Carlo wave functions. Solid line: 
analytical result. 
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Figure 4: Monte Carlo simulation of a multichannel detection experiment using the device in 
fig\^ to sample the relative phase distribution corresponding to the initial state of Fig. (a) 
Single realization of the multichannel detection: For each dephasing /3j = in/10, i = 0...9 
added to one of the input channels of the beam splitter, p+((3i) (resp. p_((3i)) particles are 
detected in the + (resp. —) output channel of the beam splitter with p+(/3i) + p-.(/3i) = p = 20. 
The obtained integers p+{f3i) (diamonds) are fitted with the function fccos 2 (0o — (3) (solid line) 
where —n/2 < <po < tt/2 is the adjustable parameter, varying from one realization to the other, 
(b) After 100 realizations of the multichannel detection (each starting with new condensates): 
histogram for the obtained values of <p . 

Figure 5: Single realization relative phase distribution at t = and at the 2 nd revival time 
t = 2tt/x for three different Monte Carlo wave functions. The parameters are as in fig^. 
From upper left to lower right the wave functions have experienced 0,3,1 and quantum jumps 
respectively. 

Figure 6: Relative phase probability distribution at t = and at the 2 nd revival time. The 
parameters are as in Fig^. Solid line: analytical prediction. Diamonds: average of 2.5 x 10 4 
Monte Carlo wave functions. 

Figure 7: Collapses and revivals of (c]jCfe) fix for a 10% asymmetry in the initial number of 
particles N a and in the condensates N a = 135.5 and Nf, = 165.5, leading to 7 a ^ 75, \i a 7^ \x\>. 
The initial total number of atoms is N = 301. The initial distribution for the difference in the 
numbers of particles n in the condensates is Gaussian with a standard deviation An = 6, and 
a non-vanishing mean value equal to 30. The other parameters are as in fig.^b. Diamonds: 
numerical result with 2.5 x 10 4 Monte Carlo wave functions. Solid line: analytical result. 

Figure 8: Collapses and revivals of (c^Cfe) fix for an initial phase state with N = 301 particles 
(solid line) and after an average over N with a Poisson distribution of parameter N = 301 
(diamonds). The effect of the average is mainly to suppress the odd revivals. The parameters 
are as in fig^b and the results are obtained from the analytical predictions. 

Figure 9: Collision fluxes (stars) and \^ (diamonds), due to one-body and three-body 
collisions respectively, calculated as in fig$, and inverse of the first revival time l/i rev = x/ 71 
(solid line) as a function of the total number of atoms. The trap frequency is Q = 2tc x 500Hz 
in (a) and Q = 2ir x 200Hz in (b). The vertical dashed line for N = 301 in (a) represents 
the conditions of fig ^b. X^ 1 ' corresponds to a lifetime due to background gas collisions of 350 
seconds. 
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